The future evolution of Unteraargletscher, a large valley glacier in the Swiss Alps, is assessed for the period 2005 to 2050 using a flowline model. Detailed measurements of surface velocity from the last decade allow us to relate ice flux to glacier thickness and width. Mass balance is calculated using a distributed temperature-index model calibrated with ice volume changes derived independently from comparison of repeated digital elevation models. The model was validated for the period 1961 to 2005 and showed good agreement between the simulated and observed evolution of surface geometry. Regional climate scenarios with seasonal resolution were used to investigate the anticipated response of Unteraargletscher to future climate changes. Three mass balance scenarios were defined, corresponding to 2.5%, 50%, and 97.5% quantiles of a statistical analysis of 16 different climate model results. We present a forecast of the future extent of Unteraargletscher in the next five decades and analyze relevant parameters with respect to the past. The model predicts a retreat of the glacier terminus of 800–1025 m by 2035, and of 1250–2300 m by 2050. The debris coverage of the glacier tongue reduces the retreat rate by a factor of three. The thinning rate increased by 50–183% by 2050 depending on the scenario applied, compared to the period 1997 to 2005.
Introduction
Glaciers are considered to be sensitive indicators of climate change (e.g. IPCC, 2001; Oerlemans and Fortuin, 1992; Haeberli, 1995). During the last century, alpine glaciers suffered major mass losses (e.g. Vincent et al., 2004). Large rearrangements in alpine systems due to glacier retreat or complete disappearance will result from the likely acceleration of climatic warming during the next decades. These changes affect the environment on local and on global scales. Water resources are expected to diminish, and significant societal impacts in peripheral regions are anticipated. Therefore, the investigation of glacier retreat and the development of methods for its prediction are important in order to be prepared for the new environmental situation.
The parameterization of glacier mass balance is of crucial importance to the investigation of glacier reaction to a warming climate. In recent years a wide range of models has been developed to calculate mass balance using climatological observations (e.g. Braithwaite, 1995; Hock, 1999; Hock and Holmgren, 2005; Schäfli et al., 2005; Pelliciotti et al., 2005; Gerbaux et al., 2005). Energy balance models directly address the physical processes at the glacier surface. Temperature-index models are based on a linear relation between positive air temperature and melt rate. Ohmura (2001) demonstrated that the physical base of temperature-index modeling is stronger than previously assumed. Melt is highly correlated with longwave heat flux, for which air temperature is a good indicator.
Several studies investigated the consequences of changing temperature and precipitation on alpine glaciers. Some works are focused on sea level change (e.g. Zuo and Oerlemans, 1997; Braithwaite and Zhang, 1999; Van de Wal and Wild, 2001). Schäfli (2005) conducted a study on the change in hydrological variables in several highly glacierized alpine catchment basins with a statistical distribution of regional climate models, but without taking the ice motion into account. The reaction of the glacier front to a shift in climate variables is related to the ice flow (Jóhannesson et al., 1989). Therefore, reasonable predictions of the response of an individual glacier with sufficient spatial and temporal accuracy are only possible using a combined ice-flow mass-balance model approach. Wallinga and van de Wal (1998) and Oerlemans et al. (1998) applied flowline models with simple parameterizations of climate change to a number of alpine glaciers. Schneeberger et al. (2001) used a 3-dimensional flow model and a glacier melt model forced by downscaled general circulation model (GCM) outputs to calculate the response of Storglaciären, Sweden.
The assessment of ice volume and glacier extent in the next decades requires the incorporation of climate change scenarios. There are several factors of uncertainty in climate projections. First of all, the socioeconomic and technological development of our civilization, which will determine the future emission of greenhouse gases, is not known. Our knowledge of the climate system and its processes and feedback mechanisms is still limited and there are various ways to describe the system by using quantitative models. In recent years, a large number of climate change models has been developed (IPCC, 2001). A probability-distribution function (PDF) can be derived by evaluating a whole set of model results with different, but realistic, boundary conditions. This has been done on global (Wigley and Raper, 2001) and on regional scales (e.g. Frei et al., 2006). A PDF provides a better quantitative estimation of the modeling uncertainties, and allows the future range of climate change to be inferred.
In this study we present a combined ice-flow mass-balance model, which is driven by regional climate scenarios for temperature and precipitation changes with seasonal resolution. The high temporal and spatial resolution used for this work allows us to perform detailed predictions of future glacier evolution based on field measurements and robust modeling approaches.
The presented study was initiated as part of an extension project of the hydroelectric power company Kraftwerke Oberhasli AG (KWO). The company manages a system of several dammed lakes in the catchment basin of the Aare river. To increase the system performance, the possibility is considered of raising the maximum lake level of Grimselsee (lake), its largest reservoir, by 23 m. This lake is situated in the valley of Unteraargletscher, 1 km downstream from the current glacier terminus. The planned raising of the lake level would flood the glacier forefield up to the current glacier snout. In order to assess whether the lake water would come into contact with the glacier after the completion of the project (in 2012), we performed numerical simulations of the future glacier evolution until 2050. To achieve this goal, we developed a flowline model applying novel approaches to ice-flow and mass-balance computation by carrying out the following steps:
the setup of an empirical ice-flow model based on flow-speed measurements,
distributed modeling of mass balance using temperature and precipitation data,
verification of the flowline model in the past by comparing the results with measured glacier surface geometries, and
incorporation of state-of-the-art climate scenarios into the mass-balance model to simulate future conditions.
Study Area and Relevant Data
With a length of approximately 13 km, Unteraargletscher is the fourth largest glacier in the Swiss Alps. The temperate valley glacier covers an elevation range from 1900 to 4000 m a.s.l. Two main tributaries, Lauteraargletscher and Finsteraargletscher, merge to form the Unteraargletscher extending about 5 km eastward of the confluence area with a mean slope of 4° (Fig. 1). Unteraargletscher is largely debris-covered and characterized by the dominant feature of an ice-cored medial moraine. The thickness of the debris cover is generally 0.1–0.2 m but increases progressively toward the glacier terminus (Sugiyama, 2003).
Figure 1
Map of the ablation zone of Unteraargletscher. The extent of debris-covered ice is shaded gray. The hachured area confines the part of the glacier on which the modeling was performed, and the dashed line marks the central flowline. Annual measurements of surface velocity are available for the cross-section profiles 1 to 6.
In the last two decades Unteraargletscher was the object of extensive field studies focusing on the geometry of the glacier and its bed (Funk et al., 1994; Bauder et al., 2003), mass balance (Bauder, 2001), basal processes (Fischer et al., 2001; Rousselot, 2006), hydrology (Schuler, 2002; Fischer et al., 2005), and dynamics (Gudmundsson, 1999; Sugiyama and Gudmundsson, 2004; Helbing, 2005).
The relevant data sets for the setup of the flowline model can be divided into five categories: (1) surface and bed topography, (2) thickness changes, (3) surface velocity measurements, (4) ablation measurements, and (5) climate data. Since 1924 annual measurements of thickness changes and surface velocities at 13 profiles have been conducted. After 1990, photogrammetric analysis of aerial photographs replaced previous field surveys, and the measurements on profiles were expanded by applying digital elevation models (DEMs) to the whole ablation area (Flotron, 1924–1998). We used the velocity measurements on Unteraargletscher at profiles 1 to 6 (Fig. 1) in the period 1989–2001 for the development of the ice-flow model. DEMs of the ablation area are available for 1990–1992, 1995–1999, 2003, and 2005 (Flotron, 1924–1998; Bauder, 2001; VAW/ETHZ, unpublished data). By comparing the elevation changes of the DEMs, volume balances were derived (Bauder, 2001). A dense network of radio-echo sounding profiles covers the ablation area of the glacier (Bauder et al., 2003). The accuracy of the bed topography is in the range of 5% of the ice thickness (Bauder, 2001). Direct ablation measurements on Unteraargletscher are available at up to 40 stakes for the period 1996–2001 (Bauder, 2001). Temperature and precipitation have been measured since 1959 at subdaily resolution and since 1989 at hourly resolution at the MeteoSchweiz weather station at Grimsel Hospiz, 1980 m a.s.l., located within 7 km from the glacier terminus.
Additional data cover the entire last century and provide possibilities for comparing the predictions of the flowline model with long data series of the past. Annual records of length change of Unteraargletscher extend back to 1871. Six DEMs of the whole glacial system were derived by evaluating old maps and aerial photographs (Bauder et al., 2007). Volume changes and rates of glacier thinning can be established for more than 125 years. A nearly constant glacier retreat has been observed since the records began. Short periods of positive mass balance (for example in the 1980s) are not revealed in the length change record. This is because Unteraargletscher has a response time of 70 to 100 years according to the volume-time scale of Jóhannesson et al. (1989) and does not reflect short-term climate variations.
Flowline Model
A site-specific one-dimensional flowline model was developed for Unteraargletscher below profile 6 (Fig. 1). The temporal evolution of the glacier surface is given by the continuity equation
where Q denotes the ice flux through a cross section, W the glacier width,
the mass-balance rate, and h the ice thickness. The x-axis follows the flowline. Equation (1) is solved explicitly with a grid spacing of Δx = 25 m along the flowline and a time step of Δt = 1 d. At each grid point and time step, Q, W, and 
are calculated as described in the following subsections.
Ice-flow Model
The ice velocity u observed on a glacier surface is composed of a component ud due to deformation of the ice mass and a component ub due to sliding over the bed (Paterson, 1994). The deformation of ice can be described by the conservation equations of mass and momentum together with a constitutive relationship accounting for ice rheology (e.g. Glen, 1955). The motion at the base including a variety of relevant processes for sliding and deformation of the bed itself is an important, but poorly understood, element of ice-flow dynamics. Under real conditions the stress field in glacier ice is very complex and can only be solved numerically (e.g. Gudmundsson, 1999). The detailed measurements of surface geometry and velocity made over a period of more than one decade in recent years offer an alternative way to calculate the ice flux along the flowline.
To incorporate two-dimensional effects into the ice-flow model, we define an empirical relationship between the ice flux through a glacier cross section and the central ice thickness by using measured data. The ice flux Q through the profiles 1–6 is calculated by integrating flow speed in cross sections.
taking y across the glacier and z perpendicular to the bed. W is given in Flotron (1924–1998), and the thickness h is obtained from radio-echo sounding and DEM. The vertical distribution of the horizontal flow speed u is calculated from measured surface velocity with an assumption of simple shearing (Paterson, 1994), The surface speed us across the profiles was measured annually from 1989 to 2001 and n = 3 is used as the exponent in Glen's flow law. We assume, that the basal sliding speed ub accounts for 50% of the surface speed based on bore-hole inclinometry on Unteraargletscher (Gudmundsson et al., 1999; Helbing, 2005). Variations in the fraction of ub would lead to a maximum sensitivity in the calculated flux of 10%, which lies in the range of uncertainty of the other variables.
The calculated fluxes through the six cross profiles for the years 1989 to 2001 enable us to parameterize ice flux with ice thickness at the center line of the glacier (Fig. 2a). For this empirical relation, a third-order fit yields good agreement for all profiles. Another empirical relation (Fig. 2b) between ice thickness and glacier width is used to obtain W in Equation (2). The values for the lowest two profiles do not match well with the fitted curve. Their influence on the computation, however, is not substantial as ice thickness and flow speed are minimal.
Figure 2
Empirical site-specific relations between central ice thickness and ice flux (a), and ice thickness and glacier width (b) for Unteraargletscher. Crosses correspond to data of individual years. Data clusters relate to the position of the profiles indicated by P1–P6 (Fig. 1). Equations of the polynomial fits are given.
Using this approach we can calculate the two-dimensional ice flux through cross sections directly from the central ice thickness (Fig. 2a) and obtain the one-dimensional flux based on the glacier width (Fig. 2b) at every point along the flowline. The ice flux needs to be prescribed as the boundary condition at the upper end of the model domain. This quantity is obtained from the empirical relationship shown in Figure 2a. Thus the upstream evolution of the glacier is taken into account in an indirect manner only. This parameterization is justified if we assume a constant mass-balance gradient over time.
Mass-balance Model
We use a distributed temperature-index model (Hock, 1999) coupled with an accumulation model for the calculation of the spatial distribution of mass balance. It elaborates on classical models using degree-day factors by varying these as a function of potential clear-sky radiation in order to account for the effects of slope, aspect, and shading. Temperature and precipitation at a nearby weather station are required field data. Surface melt rates M are computed by
where FM denotes a melt factor, rsnow/ice are radiation factors for snow and ice, and Ic is the potential clear-sky radiation. Ic is calculated with a DEM on a 25 m grid at hourly resolution.
Large parts of Unteraargletscher are debris-covered. If the coverage is continuously thicker than several centimeters, this causes a significant reduction of ablation (Kayastha et al., 2000). Stake measurements on Unteraargletscher show strong spatial variability of ice melt. To account for the effect of debris-covered ice, ablation was corrected using a constant reduction factor fdebris (Schuler et al., 2002). The melt rate M over debris-covered surfaces is calculated as
The factor fdebris is determined as the mean reduction of measured ablation at stakes on debris-covered ice in comparison to stakes on bare ice surfaces in the same elevation range. The factor fdebris is 0.5 and is supposed to remain unchanged over time. This assumption will be discussed later on.
Air temperature at every grid cell is calculated with a constant lapse rate of −6 × 10−3 °C m−1. The hourly measurements of temperature and precipitation at Grimsel Hospiz are taken to be representative of the modeling area. In this model approach, correction factors account for gauge under-catch errors and increased losses in the case of solid precipitation. Precipitation is assumed to increase linearly by 1.5% per 100 m. A threshold temperature distinguishes snow from rainfall with a linear transition range (0.5–2.5°C) of the fraction of the solid and the liquid phase (Hock, 1999). The spatial distribution of solid precipitation is corrected by taking into account the effects of snowdrift and avalanches. This is achieved by evaluating curvature and slope based on a DEM, similar to a procedure discussed in Bloeschl et al. (1991).
The applicability of temperature-index models depends on the calibration of the key parameters with field data. A site-specific relationship is established between melt rate M and positive air temperature (Hock, 2005). For the presented model there are three melt parameters to calibrate: FM, rice, and rsnow. The precipitation parameters mainly affect accumulation. The modeled area lies in the ablation zone in all years considered. Therefore, the results in this region are less sensitive to the precipitation parameters and they are set to standard values.
For the determination of an optimal set of the melt parameters we use ice volume changes derived from DEMs. Ten DEMs for the calibration period 1990–2005 are obtained from grid-based high resolution evaluation of aerial photographs. Uncertainty is estimated as ±0.3 m (Bauder, 2001). Border lines of the glacier are digitized from the aerial photograph. This step can pose some problems as the limit of the glaciated area is not always visible because of debris coverage. By comparing elevation changes between the DEMs within the glacierized area of the modeled region we calculate mean thickness changes. By assuming an ice density of 900 kg m−3 we obtain the change in water equivalent between two DEMs. The overall accuracy of the volume change evaluation is estimated as ±10% for annual time steps and ±5% for longer periods. Bauder et al. (2007) give a detailed overview about the calculation of volume changes based on DEMs and their uncertainties.
As the DEMs do not cover the entire glacierized area, the derived volume changes cannot be compared directly to the simulated mass balances. The volume change ΔVmr of the modeled region has to be corrected with the ice flux from the accumulation area counteracting the thinning of the glacier tongue due to ablation. The mean specific mass balance of the modeled region 
is therefore
The applied calibration procedure proved to be well suited for the purpose of this study, because all quantities in Equation (6) are not point based, rather they integrate overall spatial changes of the entire glacier. Our focus when calibrating was on cumulative results of simulated mass balance and corrected volume change. Thus deviations in individual years due to inaccurate determination of mass balance, ice flux, or volume change are averaged out. Good agreement was achieved between calculated mass balance and corrected volume change in the modeled region (Fig. 3). Subsequently, the results are compared to the measured melt at ablation stakes (Bauder, 2001), which serve as independent verifiers of the calibration procedure. The field data correspond well to the simulations. Mass balance is not modeled accurately at all stake locations, which can be explained by the considerable inhomogeneity of the debris coverage. However, the variability of ablation is averaged out spatially. The simulations of the mass-balance model calibrated with corrected volume changes are thus consistent with the direct ablation measurements.
Figure 3
Calibration of the model parameters with simulated mass balance (solid line) and corrected volume change (crosses) in the modeled region.
Melting in the ablation area is characterized by two different mass-balance gradients (below 2500 m a.s.l.) due to the two different ice surfaces: debris-covered, and bare (Fig. 4). The measurements at ablation stakes support this finding (Bauder, 2001). The larger scatter in the accumulation area (above ≈2900 m a.s.l.) is due to snowdrift and slope effects.
Figure 4
Distribution of mass balance with elevation: October 1998 to October 1999. The solid line represents the mean mass balance in elevation classes of 25 m and the dashed lines the root mean square-deviations of mass balance from the mean value at grid cells (dots). Below 2500 m a.s.l. there are two clearly distinguishable gradients for bare and debris-covered ice. Triangles show the measured ablation in 1998/1999.
In order to relate mass balance to the flowline with a variable geometry, we parameterize mass balance dependent on elevation. Accumulation and ablation at a given altitude is a function of climate variables and is calculated by the model described above. We assume slope and aspect at a fixed elevation to remain unchanged. Thus, the mass-balance gradient does not change, which corresponds to observations on alpine glaciers (Hoinkes, 1970). It is shifted according to the perturbations of the climate variables.
Future Climate
Climate Scenarios for the Alpine Region
The climate projections are based on simulations of the PRUDENCE project (Christensen et al., 2002). A statistical analysis from results of 16 regional climate models for Switzerland was performed by Frei (2005). These results are derived from model chains of different emission scenarios and downscaled GCM outputs. Frei (2005) computed a 95% confidence interval for the evolution of temperature and precipitation in a seasonal resolution for 2030, 2050, and 2070. Figure 5 shows the results of this study. The given values are changes relative to 1990. The temperature rise is most pronounced in the summer season and is +1.4°C until 2030 and +2.7°C until 2050, respectively. In winter, a slight increase of precipitation (+8%) until 2050 is expected, whereas a decrease of −17% until 2050 in summer is predicted (Frei, 2005).
Figure 5
Expected climate change for the northern flank of the Alps in seasonal resolution for 2030, 2050, and 2070. The bars show the 95% confidence interval, the horizontal lines the median value. Temperature (a) and precipitation (b) scenarios are relative to 1990. Figure by Frei (2005).
Mass-balance Scenarios
We define three climate change scenarios: two extreme evolutions and a median scenario, expected to be the most probable evolution. Scenario 1 is most favorable to glacier existence, Scenario 2 adopts the median of the statistical analysis (Frei, 2005), and Scenario 3 is based on the most dramatic climate change assumptions (Fig. 5 and Table 1Table 2).
Table 1
Model parameters, values and units.
Table 2
Definition of the climate scenarios according to Frei (2005) and Figure 5.
We use the model presented in an earlier section (MASS-BALANCE MODEL) to calculate future mass balances. By applying this method we do not have to rely on statistical relations between climate variables and mass balance as proposed by previous studies, but calculate accumulation and melt directly in high temporal resolution. Thus, feedback mechanisms, such as longer ablation seasons, larger portions of liquid precipitation, etc. are implicitly included. The parameters of the mass-balance model were calibrated for the period 1990 to 2005. To simulate mass balances in the next 50 years using the same parameter set, time series of temperature and precipitation with identical temporal resolution are required. In the future, the mean of climate variables is supposed to change, whereas we assume their characteristics (e.g. diurnal amplitudes) to remain the same, thus justifying the applicability of the model parameters for the next decades.
Based on each year of a 1990 to 1999 reference period, a temperature time series of hourly resolution was generated as follows: each hour value i of year j (j = 1990,…, 1999) is corrected with a normalization offset OT computed with the mean seasonal temperature of the individual year 
and the mean seasonal temperature of the reference period 
with the relation
Results
Model Validation
The model was applied to two periods in the past, for which the model output can be compared with corresponding measured glacier surface geometry. For the short period of 1990 to 2005, good agreement was found. Even for the longer period of 1961 to 2005, the simulated glacier surface matches the observed one (Fig. 6). We conclude that the combined ice-flow mass-balance model is capable of simulating glacier evolution caused by a modification of present climate input parameters over a period of a few decades. The good performance of the model between 1961 and 2005 without tuning of the parameters of both the ice-flow and the mass-balance model strongly supports the application of the flowline model for a period of comparable length in the future.
Retreat of Unteraargletscher
By forcing the model with the climate scenarios described in an earlier section (MASS-BALANCE SCENARIOS), we calculate the future evolution of Unteraargletscher until 2050. We present the results focused on 2035 and 2050. The model runs until 2035 show a retreat of the snout of Unteraargletscher between 800 m and 1025 m, depending on the scenario applied. A retreat of between 1250 m and 2300 m is predicted by 2050. Between 2035 and 2050 the rate of retreat increases substantially. Figure 7 demonstrates the glacier evolution along the central flowline. The simulated ice surfaces for the three climate scenarios corresponding to the 95% confidence interval of expected climate change all show a major retreat (Fig. 7a).
Figure 7
(a) Comparison of the model results for Scenarios 1 to 3 (2035) along the central flowline; (b) simulated glacier surface elevation for 2015, 2025, 2035, and 2050 according to Scenario 2. Shaded areas in front of the glacier snout mark the extension of the current lake (dark gray) and the projected lake (light gray).
As the climate scenarios given by Frei (2005) provide information until 2070, we conducted model runs with less reliable mass-balance assumptions due to an extremely accelerated climate change until 2070. For Scenario 1, only a slight acceleration of the retreat rate can be observed, whereas for Scenario 3 a complete disintegration of the glacier in the model domain is predicted (Fig. 8a).
Figure 8
(a) Observed and simulated length change of Unteraargletscher. The dashed line shows the trend of the measured length change extrapolated into the future. (b) Ice volumes in the modeled region and corresponding volume change rates derived from digital elevation models (DEMs) since 1880 until present and from simulation results until 2070. (c) Thinning rates derived from DEMs and simulations.
Simulated length, volume changes, and thinning rates are presented in Figure 8 with respect to observations in the past century. Ice volume in the modeled region will decrease by 41–57% by 2035 and by 60–83% by 2050 relative to 2005. However, the rate of volume change remains nearly constant (Fig. 8b). The thickness changes of Unteraargletscher were computed for six periods in the last century. Annual thinning rates in the modeled region lie in the range of −0.8 to −2.3 m a−1 with a weighted mean of −1.24 m a−1 between 1880 and 2005 (Fig. 8c). In the period 2005 to 2035, the ice thickness in the modeled region changes at an average rate of −2.7 m a−1 for Scenario 1, −3.2 m a−1 for Scenario 2, and −4 m a−1 for Scenario 3 (Table 3). The simulated thinning rates were evaluated in 5-year periods until 2050. An acceleration of thinning is predicted. Compared to the latest evaluated period, 1997 to 2005, an increase in the thinning rate of 50% by 2050 for Scenario 1, of 106% for Scenario 2, and of 183% for Scenario 3 is calculated.
Table 3
Compilation of model results for Scenarios 1 to 3 for two periods (2005 to 2035, 2035 to 2050).
We assess the mass budget of the modeled region by comparing mass gain due to ice flux through profile 6 with mass loss due to ablation. According to the applied climate scenario, the ratio of the importance of ice flux to mass balance falls from between 30 and 45% in 2005 to between 13 and 27% in 2035 and between 5 and 18% in 2050 depending on which scenario is adopted (Fig. 9). The computed mean specific mass balance 
of the modeled region of Unteraargletscher shows a progressive decrease. Because of the accelerated retreat due to climatic forcing, the importance of ice flow becomes increasingly marginal but not negligible compared to mass balance.
Figure 9
Comparison of absolute specific values (m WE) of mass fluxes in the modeled region. Ice flux through profile 6, positive, is plotted in dark gray; mean specific mass balance in the modeled region, negative, in light gray. The solid line is the ice-flux mass-balance ratio.
The computed glacier margins in the future are projected on a map assuming the thickness change on the flowline to be identical at the same elevation on the whole glacier. In Figure 10 the spatial extent of retreat is depicted. Close to the snout the glacier shrinks considerably. The model results predict that Unteraargletscher will lose 27–38% of its area downstream of profile 6 by 2035 and 40–61% by 2050 (Fig. 10 and Table 3).
Figure 10
Simulated retreat of Unteraargletscher in 2015, 2025, 2035, and 2050 for the three different climate change scenarios. Dotted transverse lines indicate the location of radio-echo sounding profiles. The dashed line in (b) follows the central flowline used for the calculations. After 2035 a proglacial lake starts to form, shown by gray patches. Shaded areas in front of the glacier snout mark the extension of the current lake (dark gray) and the projected lake (light gray). Axes are labeled in kilometers of the Swiss national coordinate system, with north at the top.
The debris coverage of Unteraargletscher delays the retreat of the glacier tongue considerably. We performed two sensitivity tests: (1) we removed the debris cover completely, or (2) we decreased fdebris [Equation (5)] at constant rates corresponding to a gradual 50% thickening of the debris cover by 2050. The results of experiment 1 show that the retreat of the glacier tongue would be about three times faster without any debris coverage compared to the reference runs. For experiment 2 we observe a deceleration of glacier retreat rate by 20%. In both experiments the perturbations of mass balances caused by different debris coverages are amplified by a positive feedback due to ice flow.
Discussion and Conclusion
Studies forecasting the future have to cope with problems of extrapolation. While our model was successfully tested for five decades in the past (Fig. 6), its performance in the future cannot be proved. This poses limitations on the time period of extrapolation into the future. The parameters of the temperature-index model were calibrated over the past 15 years by integrating the prevailing climate characteristics through this period. With major changes in the climate variables, the tuned melt parameters may be influenced over the following decades by altered circulation patterns and the redistribution of heat fluxes, which are difficult to predict (Hock, 2005). Thus, the reliability of the mass-balance predictions decrease with extreme climate change and could be subject to question by the end of the 21st century. Therefore, we confine our main conclusions to 2050.
Debris coverage is an important factor when modeling the retreat of Unteraargletscher. It is a difficult task to predict the evolution of the debris cover and its insulating effect over time spans of several decades. Long-term observations of debris thickness evolution have not yet been conducted. For Unteraargletscher we evaluated the extent of debris-covered glacier surface from topographical maps of 1880 and 1927 and from aerial photographs of the last four decades. We find that the boundaries of debris-covered ice change very slowly over time. The debris-covered surfaces shifted to higher elevations with glacier retreat, but their total area remained constant during the last century. This implies that there is some long-term balance of debris coverage although no information about debris thickness is given.
Many studies have addressed ablation on debris-covered glaciers either theoretically (e.g. Anderson, 2000) or empirically (e.g. Lundstrom et al., 1993; Kayastha et al., 2000; Sugiyama, 2003). Generally, a thickening of the debris coverage must be expected with glacier retreat. This leads to a reduction of ice melt (Lundstrom et al., 1993). However, on down-wasting debris-covered glacier tongues some processes exist that contribute to an evacuation of debris into the forefield. Lukas et al. (2005) report supraglacial debris flows exposing steep ice cliffs which melt back at an enhanced rate. Supraglacial streams also contribute to debris transport (Anderson, 2000; Lukas et al., 2005). Such meltwater channels as well as ice cliffs are presently observed on Unteraargletscher. The cliffs concentrate at the glacier terminus and along the medial moraine. The combination of all processes mentioned above inhibit well-founded forecasts of the evolution of fdebris [Equation (5)] in the next 50 years. We believe that debris thickening and thinning will balance itself out with ongoing glacier retreat. Therefore, we decided to fix the reduction factor for debris on a constant value.
The model results show retreat rates until 2035 in the range of 30 m a−1. The mean retreat rate between 1871 and the present amounts to 17.5 m a−1 (Fig. 8a). We therefore predict an increase of glacier retreat rate of approximately 70% in the next three decades compared to the past century. At first, the reaction of the glacier is rather moderate compared to anticipated increase of summer temperature given by the climate models. Homogenized data sets show a trend of ≈0.009°C a−1 between 1864 and 2004 (Begert et al., 2005). Summer temperature is expected to rise by between 0.023°C a−1 and 0.078°C a−1 by 2050, according to Frei (2005). However, length variations of a glacier are largely affected by the geometry of the ice mass and its underlying bedrock topography. The rapid retreat of the glacier snout takes place with a delay of several decades after the climate change (Fig. 8a). This lag is due to a substantial thinning of the flat glacier tongue of Unteraargletscher, which precedes a major length change.
If ice flux through profile 6 equals mass loss due to negative mass balance over the modeled region, the glacier surface remains unchanged and the position of the snout is stationary, as long as no modifications of material and basal properties occur. At present, the ratio ice-flux to mass-balance in the modeled region lies in the range of 40%, indicating a pronounced imbalance of the mass budget. Under the influence of climate variations this ratio can change rapidly, as it is directly influenced by the glacier mass balance (Fig. 9). The boundary condition at the upper end of the model, i.e., the starting point of the flowline, is the ice flux. Assuming a changing mass-balance gradient, this may pose some problems. With large changes in climate variables, the mass-balance gradient is likely to be altered and, hence, some of our assumptions may no longer be justified. However, the importance of ice flux compared to mass balance decreases with time and approaches 10% by 2050 (Fig. 9). Thus, the flux at the upper model boundary loses its influence on the final results in a gradual way. Additionally, the test runs (Fig. 6) indicate that the parameterization of the boundary condition at the upper end of the model domain reproduces reasonable results, also for conditions that are not steady state.
Ice flux should not be neglected in glacier retreat studies. Although the influence of ice melt exceeds the compensating effect of ice flow by far in the case of Unteraargletscher in the 21st century, significant underestimation of glacier extent occurs when assuming solely non-dynamic down-wasting of the ice. Test runs excluding ice flow show that the mean thinning rate is overestimated by 10% to 25% by 2050 with larger errors for Scenario 1 than for Scenario 3.
Ice volume is an important parameter for the assessment of water resources. The total ice volume in the modeled region currently amounts to 0.75 km3, corresponding to a decrease of 60% since 1880 (Fig. 8b). For all climate scenarios ice volume in the modeled region shrinks steadily, but not with an accelerated trend. The rate of volume change even shows an upward trend in future in spite of continuing climate warming. Due to the progressive decrease of glacier area, the meltwater production will be reduced as well. This is a key observation for the management of water resources in catchment areas dominated by glaciers. The future meltwater discharge from highly glacierized drainage basins is expected to increase, favored by additional ablation in a first stage after the shift of climate variables. As ice volume and glacier area shrink, meltwater discharge drops below the preceding level. Thus, the short-term (daily to annual scales) and long-term (decadal scales) storage capacities of glaciers are lost. This will almost certainly lead to a reduction of water resources during dry and hot summers (Hock et al., 2004).
This work shows the impact of expected climate change on the evolution of a typical alpine glacier. Detailed simulations of future glacier mass balance were performed using a probability-distribution function of different regional climate scenarios for temperature and precipitation change in seasonal resolution. By combining the calculations of ice flow and mass balance, predictions of the glacier extent in the next 50 years are possible.
The unique database existing for Unteraargletscher allows an empirical ice-flow model to be set up based on velocity and ice thickness measurements and optimal calibration and verification facilities for the mass-balance model. According to the results of this study, the projected lake will never reach the glacier terminus. In the next decades, the glacier snout will retreat substantially and considerable growth of the proglacial area is expected. The retreat and thinning rates will double by the middle of the 21st century, leading to a dramatic decrease of ice volume.
Acknowledgments
The study was commissioned by Kraftwerke Oberhasli AG (KWO). Flotron AG made available surface velocities at different profiles and evaluated digital elevation models. We acknowledge C. Frei for the climate scenario data. R. Hock provided the code of her model used to calculate the mass balances. We thank H. Blatter for valuable comments on an earlier version of the manuscript. Suggestions of the scientific editor S. Anderson and reviews of B. Hubbard and an anonymous reviewer helped to improve the paper.
